function [Q,F,S,empty] = drawQ(restr,irs,phi,opt)
% Function attempts to draw Q from the space of orthonormal matrices 
% satisfying the equality and sign restrictions using Algorithm 1 (Step 2).
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options

eqRestr = restr.eqRestr;
signRestr = restr.signRestr;
f = restr.f;
s = restr.s;
vma = irs.vmaGK;
lrcir = irs.lrcirGK;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
B = phi.B;
p = opt.p;
const = opt.const;
signNorm = opt.signNorm;

N = size(Sigmatr,1); % Number of variables in VAR
mZ = size(eqRestr,1); % Number of equality restrictions
mS = size(signRestr,1); % Number of sign restrictions

% Construct matrix representing equality restrictions.
% F(phi) = [F_1(phi); ...; F_N(phi)], where F_i(phi) represent equality 
% restrictions on the ith column of Q, so F_i(phi)*q_i = 0.

if mZ == 0 % If no equality restrictions

    F = [];

else % If equality restrictions

    F = zeros(mZ,N);

    for ii = 1:mZ

        if eqRestr(ii,4) == 1 % Restriction on A0

            F(ii,:) = Sigmatrinv(:,eqRestr(ii,2))';

        elseif eqRestr(ii,4) == 2 % Restriction on A0^(-1)

            F(ii,:) = Sigmatr(eqRestr(ii,1),:);

        elseif eqRestr(ii,4) == 3 % Restriction on A_l
            
            B = reshape(B,[N*p+const,N]);
            Bl = B(const+(eqRestr(ii,3)-1)*N+1:const+(eqRestr(ii,3))*N,...
                eqRestr(ii,2));
            F(ii,:) = (Sigmatrinv*Bl)';
            
        elseif eqRestr(ii,4) == 4 % Restriction on LRCIR

            F(ii,:) = lrcir(eqRestr(ii,1),:);

        end

    end

end

% Order rows of F in terms of the column of Q restricted.

Ftilde = [];

if isempty(F)
    
    Ftilde = F;
    
else

    for jj = 1:N % For each column of Q

       % Find rows of F restricting jth column of Q.
       Fj = F(eqRestr(:,1)==jj & (eqRestr(:,4)==1 | eqRestr(:,4)==3)...
           | eqRestr(:,2)==jj & (eqRestr(:,4)==2 | eqRestr(:,4)==4),:);
       Ftilde = [Ftilde; Fj];

    end
    
end

% Construct matrix representing inequality restrictions (excluding 
% normalisations).
% S(phi) = [S_1(phi); ...; S_N(phi)], where S_i(phi) represent inequality 
% restrictions on the ith column of Q, so S_i(phi)*q_i >= 0.

S = zeros(mS,N);

for ii = 1:mS
    
    if signRestr(ii,5) == 1 % Sign restriction on IRF
        
        S(ii,:) = vma(signRestr(ii,1),:,signRestr(ii,3)+1)*signRestr(ii,4);
        
    elseif signRestr(ii,5) == 2 % Sign restriction on A0
        
        S(ii,:) = Sigmatrinv(:,signRestr(ii,2))'*signRestr(ii,4);
        
    end
    
end

% Order rows of S in terms of the column of Q restricted.

Stilde = [];

if isempty(S) 
    
    Stilde = S;
    
else
    
    for jj = 1:N % For each column of Q

       % Find rows of S restricting jth column of Q.
       Sj = S((signRestr(:,1)==jj & signRestr(:,5)==2) | ...
           (signRestr(:,2)==jj & signRestr(:,5)==1),:);
       Stilde = [Stilde; Sj];

    end
    
end

iter = 0;
empty = 1;

while iter <= opt.maxDraw && empty == 1

    iter = iter + 1;

    % Draw Q from space of orthonormal matrices satisfying equality 
    % restrictions.

    % Generate vector of standard normal random variables.
    z = randn(N,1); 

    qtilde = zeros(N);
    
    if isempty(F)

        qtilde(:,1) = z;
        
    else
    
        % Find rows of F restricting first column of Q (F1).
        F1 = Ftilde(1:f(1),:);

        % Compute residual from linear projection of z on F1'.
        [q,r] = qr(F1',0);
        qtilde(:,1) = z-F1'*(r \ (q'*z));

    end

    for jj = 2:N % For each column of Q

        if isempty(F)
            
            regMat = qtilde(:,1:jj-1);
            
        else
        
           % Find rows of F restricting jth column of Q (Fj).
           Fj = Ftilde((sum(f(1:jj-1))+1):sum(f(1:jj)),:);
           regMat = [Fj' qtilde(:,1:jj-1)];
        
        end

       % Generate vector of independent standard normal random variables.
       z = randn(N,1); 

       % Compute residual from linear projection of z on regMat.
       [q,r] = qr(regMat,0);
       qtilde(:,jj) = z - regMat*(r \ (q'*z));

    end

    Q = zeros(N);

    for jj = 1:N

        if signNorm == 0

            % Normalise diagonal elements of A0 to be positive.

            if Sigmatrinv(:,jj)'*qtilde(:,jj) == 0

                % Randomly set sign of q_j to -1 or 1.
                Q(:,jj) = (1-2*(rand > 0.5))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

            else

                Q(:,jj) = sign(Sigmatrinv(:,jj)'*qtilde(:,jj))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

            end 

        else % Normalise diagonal elements of A0^(-1) to be positive

           if Sigmatr(jj,:)*qtilde(:,jj) == 0

                % Randomly set sign of q_j to -1 or 1.
                Q(:,jj) = (1-2*(rand > 0.5))*...
                    (qtilde(:,jj)./norm(qtilde(:,jj)));

           else

               Q(:,jj) = sign(Sigmatr(jj,:)*qtilde(:,jj))*...
                   (qtilde(:,jj)./norm(qtilde(:,jj)));

           end

        end

    end
    
    % Check whether proposed draw satisfies sign restrictions (if any).

    if mS > 0  % If sign restrictions

        signCheck = zeros(mS,1);
        restrCount = 0;

        for ii = 1:N
            
            for jj = 1:s(ii)
                
                restrCount = restrCount + 1;

                signCheck(restrCount) = S(restrCount,:)*Q(:,ii);
                
            end
            
        end

        % If IRs satisfy sign restrictions, terminate while loop and 
        % record that set is non-empty.
        empty = (min(signCheck) < 0);
        
    else % If no sign restrictions
        
        empty = 0;
        
    end
    
end

F = Ftilde; % Replace F with reordered F
S = Stilde; % Replace S with reordered S

end